#include "SimulatorParameters.h"
#include <sstream>

namespace PRS{
	SimulatorParameters::SimulatorParameters(){
		TimeStepCounter = 0;
		exportIter = 0;
		accSimTime = .0;
		stop_simulation = true;
		useHOApp = false;
		restart = false;
		vtk_step = 0;

		PVI_increment = 0.01;	// every 5% of total simulation time a new VTK file will be printed
		PVI_accumulated = .0;	// summation of all PVI_increments

		allowPrintingVTK = false;
		pctype = PCNONE;
		firstVTKupdate = true;
		EBFV1_pressureSolver_scheme = false;
	}

	SimulatorParameters::~SimulatorParameters(){
		MIter_RockProperties mit = mapRockProp.begin();
		for (; mit != mapRockProp.end(); mit++){
			double *K = mit->second->K;
			delete[] K; K = 0;
		}
		mapRockProp.clear();
		mapBC.clear();
	}

	void SimulatorParameters::getNodesWithOutWells(TMesh* pTMesh_data){
		// get vertices without wells
		std::set<VertexData*> verticesWithOutWells;
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();
		for (int i=0; i<ndom; i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				if (!isInjectionWell(vdata->flag) && !isProductionWell(vdata->flag)){
					verticesWithOutWells.insert(vdata);
				}
			}
		}
		cout << "Num vertices without wells: " << verticesWithOutWells.size() << endl;
		int row = 0;
		nrows = (int)verticesWithOutWells.size();
		pTMesh_data->vertices_WithOutWells.allocateMemory(nrows,1);
		for(std::set<VertexData*>::iterator iter = verticesWithOutWells.begin(); iter!= verticesWithOutWells.end(); iter++){
			pTMesh_data->vertices_WithOutWells.setValue(row++,0,*iter);
		}
		verticesWithOutWells.clear();

	}

	void SimulatorParameters::checkIfRankHasProductionWell(TMesh* pTMesh_data){
		// check if rank has a production well
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();
		for (int i=0; i<ndom; i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				if (isProductionWell(vdata->flag)){
					setRankHasProductionWell();
					break;
				}
			}
		}

	}

	void SimulatorParameters::checkPermeabilityTensor(){
		const double* K = getPermeability(3300);//BACALHO
		// Kxx*Kyy >= Kxy*Kyx
		if (K[0]*K[3] < K[1]*K[2]){
			throw Exception(__LINE__,__FILE__,"Permeability tensor must obey the following relation: K_xx*K_yy >= K_xy*K_yx\n");
		}
		K_Isotropic = (K[1]*K[2] == .0)?true:false;
	}

	const double* SimulatorParameters::getPermeability(const int &dom){
		MIter_RockProperties mit = mapRockProp.find(dom);
		if (mit != mapRockProp.end()) return mit->second->K;
		cout << "Warning: no permeability tensor associated to domain " << dom << ".\n";
		cout << __FILE__ << "\t at line " << __LINE__ << endl;
		return 0;
	}

	void SimulatorParameters::getDomains(){
		for (MIter_RockProperties mit=mapRockProp.begin(); mit!=mapRockProp.end(); mit++) setOfDomains.insert(mit->first);
	}
		
	double SimulatorParameters::getBC_Value(const int &flag){
		MapFlagIter mIter = mapBC.find(flag);
		if (mIter != mapBC.end()) return mIter->second->val;
		cout << "Warning: getBC_Value() return null value\n";
		cout << __FILE__ << "\t at line " << __LINE__ << endl;
		return 0;
	}

	bool SimulatorParameters::isNodeFree(const int &flag){
		MapFlagIter mIter = mapBC.find(flag);
		if (mIter == mapBC.end()) return true;
		else{
			BdryConditionType *bct = mIter->second;
			return ( !bct->type.compare("dirichlet") )?false:true;
		}
	}

	bool SimulatorParameters::finishSimulation(){
		if (!stop_simulation || (getAccumulatedSimulationTime() >= getSimTime()) ){
			if (!0){
				std::cout << "#################################\n";
				std::cout << "Simulation " << (double)100.0*getAccumulatedSimulationTime()/getSimTime() << "% concluded.\n";
				std::cout << "End of Simulation\n";
				std::cout << "#################################\n\n\n";
			}
			return 1;
		}
		else{
			if (!0) std::cout << "Simulation " << (double)100.0*getAccumulatedSimulationTime()/getSimTime() << "% concluded.\n";
			return 0;
		}
	}

	// physical parameters
	double SimulatorParameters::getPorosity(const int &dom){
		MIter_RockProperties mit = mapRockProp.find(dom);
		if (mit != mapRockProp.end()) return mit->second->porosity;

		cout << "Warning: porosity() return null value.";
		cout << "Domain "<< dom << " not found.\n";
		cout << __FILE__ << "\t at line " << __LINE__ << endl;

		return 0;
	}

	void SimulatorParameters::setInitialOilVolume(TMesh *pTMesh_data){
		VertexData *vdata;
		double TPV = .0;	// Total Porous Volume
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();
		for (int i=0; i<ndom; i++){
			double TPV_tmp = .0;
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				TPV_tmp += vdata->volume/(vdata->numRC + 1.0);
			}
			TPV_tmp *= getPorosity(vdata->flag);
			TPV += TPV_tmp;
		}
		double Sw = 1.0 - Sw_initial();
		//_IOV = P_getSumDbl(TPV)*Sw;
		//todo:: parallel
		_IOV = TPV*Sw;
	}

	void SimulatorParameters::correctTimeStep(double &timeStep){
		if (firstVTKupdate){
			updatePrintOutVTKFrequency();
		}
		double timeFrequency = getPrintOutVTKFrequency();
		double accST = timeStep + getAccumulatedSimulationTime();
		if ( accST > timeFrequency ){
			timeStep = timeStep - (accST - timeFrequency);
			accST = timeFrequency;
			allowPrintingVTK = true;
		}
	}

	void SimulatorParameters::printOutVTK(TMesh* pTMesh_data, pFunc_PrintVTK printVTK){
		//allowPrintingVTK=true;
		#ifdef _SEEKFORBUGS_
			allowPrintingVTK = true;
		#endif

		#ifdef CRUMPTON_EXAMPLE
			allowPrintingVTK = true;
		#endif

		if (allowPrintingVTK){
			static int theStep = getStepOutputFile();
			char fname[256];
			sprintf(fname,"%s__%d-of-%d__step-%d.vtk",expofName.c_str(),0,1,++theStep);
			PetscPrintf(PETSC_COMM_WORLD,"\nVTK output: %s\n",fname);
			printVTK(pTMesh_data,fname);
			updatePrintOutVTKFrequency();
			allowPrintingVTK = false;
		}
	}

	void SimulatorParameters::updatePrintOutVTKFrequency(){
		firstVTKupdate = false;

		// a new VTK file will be printed at each 0.05 PVI (5% of total simulation time)
		PVI_accumulated += getPVIincrement();

		// when next vtk file must be print out
		vtk_time_frequency = PVI_accumulated*getSimTime();
	}
}
